Nonlinear waves in disordered chains: probing the limits of chaos and spreading 
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We probe the limits of nonlinear wave spreading in disordered chains which are known to localize 
linear waves. We particularly extend recent studies on the regimes of strong and weak chaos during 
subdiffusive spreading of wave packets [EPL 91, 30001 (2010)] and consider strong disorder, which 
favors Anderson localization. We probe the limit of infinite disorder strength and study Frohlich- 
Spencer- Wayne models. We find that the assumption of chaotic wave packet dynamics and its impact 
on spreading is in accord with all studied cases. Spreading appears to be asymptotic, without any 
observable slowing down. We also consider chains with spatially inhomogeneous nonlinearity which 
give further support to our findings and conclusions. 

PACS numbers: 05.45.-a, 05.60.Cd, 63.20.Pw 



I. INTRODUCTION 

Anderson localization [l[ was discovered 50 years ago 
in disordered crystals as an accumulation of single parti- 
cle electronic wavefunctions and can be interpreted as 
an interference effect between multiple scatterings of 
the electron by random defects of the potential. As 
a consequence eigenstates are no longer spatially ex- 
tended, but exponentially localized. Anderson localiza- 
tion is a universal phenomenon of wave physics, unre- 
stricted to quantum mechanics. Experimental observa- 
tions were made in noninteracting Bose-Einstcin con- 
densates (BEC) expanding in random optical potentials 
@j 9i light propagation in spatially random nonlinear 
optical media [U, @ , and in microwave cavities filled with 
randomly distributed scatterers @. Anderson localiza- 
tion is a linear wave effect, i.e. it is well-established for 
wave equations which are linear in the wave amplitude. 
However, in many cases one is confronted with a nonlin- 
ear response of the wave-carrying medium; for instance, 
large light intensities induce a nonlinear response of the 
optical medium. Electron-electron and electron-phonon 
interactions also result in substantial deviations from An- 
derson localization in solids. In experiments of Bose- 
Einstcin condensates the interatomic interactions are al- 
ways present, although they can be diminished by cither 
decreasing atomic densities or by exploiting magnetically 
tunable Feshbach resonances. 

From a mathematical perspective, a linear wave equa- 
tion is integrable, with each normal mode evolving inde- 
pendently in time. A localized wave packet in the pres- 
ence of Anderson localization will therefore stay localized 
as time evolves. Nonlinearity will usually destroy the in- 
tegrability of a system and induce mode-mode interac- 
tions. It was observed numerically that wave packets in 
such nonlinear disordered wave equations delocalize in 
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time without respecting Anderson localization limits 0- 
[Tlj . Thus, there are several intriguing questions which 
have attracted much attention during the last few years: 
(i) will Anderson localization be destroyed by arbitrary 
small strength of nonlinearity or is there a threshold be- 
low which the localization is restored? (ii) will wave 
packet spreading, if observed, last forever or will it stop 
at certain (though probably very large) time? (iii) is the 
shape of the initial wave packet crucial for the details of 
spreading? We will mainly address question (ii) here. 

Johansson ct al. [l2j conjectured that spreading must 
eventually stop and dynamics will become close to reg- 
ular, assuming that in these limits the Kolmogorov- 
Arnold-Moscr (KAM) theorem is applicable, i.e. that 
for small wave density regular nonergodic phase space 
structures predominate and the dynamics develops along 
KAM tori. Other attempts consist in a numerical scal- 
ing analysis, in order to predict and extend results 
beyond computational ability (l3| . Analytical studies 
perform perturbation theory to higher order by treat- 
ing the strength of nonlinearity as a small parameter 
|14| . conflicting with the explosive growth of secular 
terms in higher orders of perturbation theory. This 
theory states that for the disordered discrete nonlinear 
Schrodinger model with nonlinearity strength exceed- 
ing a finite threshold, any initial localized wave packet 
cannot fully spread to zero amplitudes at infinite time. 
In this part of the excitation is selftrapped as 

a result of nonlinearity induced frequency shifts, which 
tune a localized excitation out of resonance with its sur- 
rounding non-excited linear modes. However, even in 
the case of strong nonlinearity, subdiffusion of the non- 
sclftrapped part is observed [{§. When strong nonlin- 
caritics are avoided numerical studies showed a rather 
universal asymptotic subdiffusive spreading of initial sin- 
gle site excitations [!, [H UH , which is characterized by a 
growth of the second moment of the wave packet as t a , 
with q < 1 [lH . The selftrapping theorem [l6[ holds 
irrespective of the strength of disorder, therefore it is re- 
flecting the properties of a strongly nonlinear lattice wave 
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equation, rather than peculiarities of waves propagating 
in disordered media. Also the selftrapping theorem is 
crucially depending on the presence of at least two inte- 
grals of motion, and fails for most nonlinear wave equa- 
tions with only one integral of motion. 

In Q the observed wave packet spreading was assumed 
to be due to an incoherent excitation of the wave packet 
exterior, induced by the chaotic dynamics of the wave 
packet interior. The number of resonant modes in the 
packet was estimated by conside ring from quadruplet and 
triplet mode-mode interactions [17] . A generalization to 
higher dimensions D and different nonlincarity powers 
were performed. This led to a quantitative prediction for 
the subdiffusive wave packet spreading characteristics a 



|9| . Its validity was confirmed numerically in [S , 15 , [18| ■ 
Recently, it has been predicted theoretically [19|, [20| and 
verified numerically [19( that a potentially long-lasting 
strong chaos regime induces faster (though still subdiffu- 
sive) spreading, which is followed by the asymptotic and 
slower weak chaos subdiffusive spreading. Notably, pub- 
lished numerical data did not report on a further slowing 
down of spreading, when starting from the weak chaos 
regime. 

In this paper we present results of extensive numerical 
studies of wave packet spreading in various models of dis- 
ordered nonlinear one-dimensional lattices. In particular, 
we consider different initial excitations and scan the pa- 
rameter space of disorder strength and nonlinearity over 
a wide region. The main aim is to test the applicability 
of previously derived spreading laws, and to search for 
indications of a continuation of the weak chaos spread- 
ing, or for indications of a slowing down, as conjectured 
by others. 



II. MODELS 



Discrete Nonlinear Schrodinger and 
Klein-Gordon chains 



In our study we consider various one-dimensional lat- 
tice models. The first one is the disordered discrete 
nonlinear Schrodinger equation (DNLS) described by the 
Hamiltonian function 

n D = 5>|^l 2 + f IV^I 4 - (V-rnVv* + ^t+M (i) 



in which ipi are complex variables, I are the lattice site in- 
dices and p > is the nonlinearity strength. The random 
on-site energies e; are chosen uniformly from the interval 
[ — T"> T'] i with W denoting the disorder strength. The 
equations of motion are generated by tpi = &Hd / d{iip^): 



(2) 



This above set of equations conserves both the energy of 
Eq. Q, and the norm S = J2i IV'il 2 - 



The second model we consider is the quartic Klein- 
Gordon (KG) lattice, given as 
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where ui and pi respectively are the generalized coordi- 
nates and momenta on site Z, and et arc chosen uniformly 
from the interval [^,f]- The equations of motion are 
■ill = —dTix / dui and yield 
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-eiui 



— (ui+i + - 2u{) . (4) 



This set of equations only conserves the energy of Eq. ^ . 
The scalar measure of energy resulting from Eq. wc 
shall henceforth label as H . This scalar value H > 
serves as a control parameter of nonlinearity, similar to 
P for the DNLS case. 

For p = and xpi = Ai exp(-iAt), Eq. reduces to 
the linear eigenvalue problem 



XAi = eiAi - - A 



i+i 



(5) 



The normalized eigenvectors A v ^ (Y^i^-tl = 1) are the 
corresponding normal modes (NMs) , and the eigenvalues 
A„ are the frequencies of these NMs. The width of the 
cigenfrequency spectrum A„ in Eq. ([S]) is Ajy = W + 4 
with \ y e [-2-^,2+^]. The coefficient 1/(2W) in 
Eq. §3§ was chosen so that the linear parts of the Hamilto- 
nians, Eqs. (|H3p would correspond to the same eigenvalue 
problem. In the limit H — > (in practice by neglect- 
ing the nonlinear term uf/4) the KG model of Eq. ([3]) 
- with ui = Ai exp(zwt) - is reduced to the same linear 
eigenvalue problem of Eq. ([S]), under the substitutions 
A = Wuj 2 - W - 2 and e ; = W{e t - 1). The width of 
the squared frequency uj 2 spectrum is Ax = 1 + ^ with 
loI e [|, | + . Note that A D = WA K . As in the case 
of DNLS, W determines the disorder strength. 

The asymptotic spatial decay of an eigenvector is given 
by A Ut i ~ e -i /^ A "' where £(Aj/) is the localization length. 
In the case of weak disorder, W —> 0, the localiza- 
tion length is approximated [13, HU as £(A„) < £(0) ~ 
100/W 2 . The NM participation number p v = 1/ J2i A t,i 
characterizes the spatial extend of the NM. An average 
measure of this extent is the localization volume V , which 
is of the order of 3.3£(0) w 330/VF 2 for weak disorder and 
unity in the limit of strong disorder, W — > oo [l7j. The 
average spacing of eigenvalues of NMs within the range 
of a localization volume is then d w A/V, with A being 
the spectrum width. The two frequency scales d < A 
determine the packet evolution details in the presence of 
nonlincarity. 

In order to write the equations of motion of Hamil- 
tonian (Q} in the normal mode space of the system we 
insert V; = Ei/^J^ m @> with |</>„| 2 denoting the 
time-dependent amplitude of the i^th NM. Then, using 
Eq. §5§ and the orthogonality of NMs the equations of 
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motion (|2|) read 

tyv = K<t>v + fi ^2 ^vu^vaK^^va^vs ( 6 ) 

with the overlap integral 

I 

The frequency shift of a single site oscillator induced by 
the nonlinearity is Si = fi\ipi\ 2 for the DNLS model. The 
squared frequency shift of a single site oscillator induced 
by the nonlinearity for the KG system is Si = (3£';)/(2e;), 
with Ei being the energy of the oscillator. Since all NMs 
are exponentially localized in space, each effectively cou- 
ples to a finite number of neighbor modes. The non- 
linear interactions are thus of finite range; however, the 
strength of this coupling is proportional to the norm (en- 
ergy) density for the DNLS (KG) model. If the packet 
spreads far enough, we can generally define two norm 
(energy) densities: one in real space, m = \ipi\ 2 (Ei) and 
the other in NM space, n u = \(f>u\ 2 {E v ). Averaging over 
realizations, no strong difference is seen between the two, 
and therefore, we treat them generally as some character- 
istic norm (n) or energy (E) density. The frequency shift 
due to nonlinearity is then So ~ fin for the DNLS model, 
while the squared frequency shift is 8k ~ 3E/2 for the 
KG lattice. The basic characteristics of both models are 
summarized in Tabic HI 

We order the NMs in space by increasing value of the 
center-of-norm coordinate X v = ^ I A 2 l [t| [ID, [l8llT9j. 
For DNLS we follow normalized norm density distribu- 
tions z v = |0iy| 2 /X^ I^mI 2 ' while for KG we follow nor- 
malized energy density distributions z v = E v j E^ 

with E v = A 2 /2 + lu 2 A 2 /2, where A v is the ampli- 
tude of the vth. NM and uj 2 its squared frequency. Wc 
measure the second moment "^2 = — v) 2 z v (with 

v = ^2, v vz v \ which quantifies the wave packet's spread- 
ing width; the participation number P = l/^2 v Z„, i- e - 
the number of the strongest excited modes in z v ; and 
the compactness index £ = P 2 /m2, which quantifies the 
inhomogeneity of a wave packet. Thcrmalized distribu- 
tions have £ ~ 3, while ( « 3 indicates very inhomo- 
geneous packets, e. g. sparse (with many holes) or par- 
tially sclftrapped ones (see [15( for more details). In ad- 
dition, following Anderson's definition of localization [l[ , 
we measure the fraction Sv (Hy) of the wave packet 
norm (energy) in a localization volume V around the ini- 
tially excited state in real space. For a localized state 
this fraction asymptotically tends to a constant nonzero 
value, while it goes to zero in the case of derealization. 



B. Frohlich-Spencer- Wayne chain 

In the limit of strong disorder (W —> 00) the DNLS and 
KG models suffer from increasing computational times 



needed to observe any nontrivial dynamics. This is be- 
cause the eigenvectors tend to single site profiles, i.e. the 
overlap integrals become very small. Frohlich, Spencer 
and Wayne (FSW) suggested considering a modified 
Hamiltonian, which operates directly in normal mode 
space for the strong disorder limit, but considers artifi- 
cial rcscalcd anharmonic interactions between neighbor- 
ing NMs in order to rescale time (22j: 

where the NMs are equivalent to the single site oscilla- 
tors. The NM eigenvalues e„ are considered to be uncor- 
rected, also for nearest neighbors. This is different from 
the DNLS and KG models. Also the FSW chain has only 
pair interactions between NMs (sites) . Note also that the 
nonlinear part of the FSW Hamiltonian is invariant un- 
der any shift u v — > u v + a, as opposed to the KG model. 

C. Models with spatially inhomogeneous 
nonlinearity 

We also consider two variants of DNLS and KG mod- 
els with spatially inhomogeneous nonlinearity. The first 
type of lattices is composed of linear coupled oscillators 
except for a central region of length L where nonlineari- 
ties are present. We refer to this type as the LNL (Linear- 
Nonlinear-Lincar) model. The second type is called NLN 
(Nonlinear-Lincar-Nonlinear) and is the exact counter- 
part of the previous one, since the linear part of the lat- 
tice is located at the central L sites. 



III. WAVE PACKET EVOLUTION 
A. Theoretical predictions 

1. DNLS and KG 

The evolution of wave packets in nonlinear disordered 
chains can be expected to be selftrapped for strong non- 
linearities, or show no selftrapping for weaker nonlinear- 
ities. The existence of the selftrapping regime was theo- 
retically predicted for the DNLS model in [l6| (see also 
fla for more details). According to the theorem stated 
in [111, for large enough nonlinearities (Sd > A-d) single 
site excitations cannot uniformly spread over the entire 
lattice. Consequently, a part of the wave packet will re- 
main localized, although the theorem does not prove that 
the location of this inhomogeneity is constant in time. 

If the nonlinear shift S moves the frequencies of some 
of the initially excited oscillators out of the linear spec- 
trum, it tunes them out of resonance and part of the wave 
packet will be sclftrapped. In our study we consider ini- 
tial "block" wave packets, where L central oscillators of 
the lattice are excited having the same norm (energy). 
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DNLS 


KG 


On-site energies 


2 2 


"1 3" 

2 2 


Spectrum 


W W~ 


2 [13 4 " 
W " e 2'2 + W 


Spectrum width A 


A D = W + 4 




Localization volume V < 


W ->■ 
W — > oo 


330 

~ r 


V ~ 1 


Average spacing d < 


• 

W -> 
W — > oo 








djc ~ const. 


Nonlinear energy shift S 


5 D ~ /3n 





TABLE I: Characteristic quantities of the DNLS |T]) and the KG (|3} models. The dependence on the strength of disorder W 
of both the localization volume V and of the average spacing d of NMs' eigenvalues within the range of V, is given for the 
limiting cases of weak W — > and strong disorder W — ► oo. Note that n (E) represents a general characteristic norm (energy) 
of wave packets of the DNLS (KG) model. 



Since we consider many random disorder realizations (of 
the order of few hundreds) we expect that, on average, 
the linear frequencies of the initially excited lattice sites e; 
(?;) cover the whole range of permitted values [— ^r] 
([|, II). Thus, some of these frequencies are tuned out 
of resonance if 5d > 2 (5d > 1/W). These conditions 
for the possible appearance of selftrapping are less strict 
than the theoretically defined ones [15|, Ha] and are, in 
general, in good agreement with numerical simulations. 
In particular, the selftrapping regime was numerically ob- 
served for single-site excitations 0, [HI [HI and for ex- 
tended excitations [l9|, both for the DNLS and the KG 
models, despite the fact that the KG system conserves 
only the total energy E, and the selftrapping theorem 
can not be applied there. 

When selftrapping is avoided for Sd < 2 (5k < 1/W), 
two different spreading regimes were predicted, having 
different dynamical characteristics: an asymptotic weak 
chaos regime, and a potential intermediate strong chaos 
one [2(j ■ Numerical verifications of the existence of these 
two regimes were presented in 0, [l9[. In the weak 
chaos regime, for L > V and 6 < d most of the NMs 
are weakly interacting with each other. Then the sub- 
diffusive spreading of the wave packet is characterized 
by TO2 ~ i 1 / 3 . If the nonlinearity is weak enough to 
avoid selftrapping, yet strong enough to ensure S > d, 
the strong chaos regime is realized. Wave packets in this 
regime initially spread faster than in the case of weak 
chaos, with 11J2 ~ t 1 / 2 . Since the norm density drops 
with further spreading, S is dropping in time as well, and 
eventually the wave packet enters the weak chaos regime, 
where its evolution is characterized by slower spreading 
with 777,2 ~ i 1 / 3 . The wave packet evolution in both the 
weak and the strong chaos spreading regimes is also ex- 



pected to be characterized by an increase of the partici- 
pation number as P ~ i"/ 2 when 7772 ~ t a . 

Let as now discuss the spreading of wave packets when 
L < V. The packet will initially spread over the localiza- 
tion volume V during a time interval Ti n ~ 2n/d, even in 
the absence of nonlinearities [2(| ■ The initial average 
norm (energy) density rii n (Ei n ) of the wave packet is 
then lowered to 77(T in ) ps n in L/V (E(T in ) rj E in L/V). 
The further spreading of the wave packet in the presence 
of nonlinearities is then determined by these reduced den- 
sities. Note that for single-site excitations (L = 1), the 
strong chaos regime therefore completely disappears and 
the wave packet evolves either in the weak chaos or in 
the selftrapping regimes [1, B EE HH • 



2. FSW chain 

There are no existing theoretical predictions for wave 
packet dynamics in the FSW chain. The FSW case can 
be considered as a strong disorder limit of the KG model, 
emulating the dynamics of the latter in NM space. How- 
ever, in the KG and DNLS case, triplet interactions be- 
tween NMs are present and necessary in order to allow 
for finite (though small) resonance probabilities for small 
(but finite) energy and norm densities [13, H3| • Pair in- 
teractions will cease to produce NM resonances for suffi- 
ciently small densities due to level repulsion within one 
localization volume [U, E3, Hil . The FSW chain keeps 
only pair interactions. At the same time, level repulsion 
between neighboring (interacting) NMs are absent in the 
FSW case. Therefore, a theoretical analysis similar to 
the DNLS and KG case [2(| appears to be possible. Its 
details will be considered in a future publication. The 
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width of the linear spectrum Apsw = 1- The average 
spacing of nearest neighbor eigenvalues dpsw ~ ^fsw- 
Therefore, we expect that only the asymptotic regime of 
weak chaos, and the regime of selftrapping can be ex- 
pected. 



3. LNL and NLN chains 

The LNL chain can be expected to start in a chaotic 
wave packet spreading regime as long as the wave packet 
is confined mainly to the finite size N (nonlinear) part. 
However, the more the wave packet spreads, the more 
it extends into the infinitely extended L (linear) parts. 
Resonances and chaos are therefore confined to the finite 
N part. Since distant NMs in the L part are exponen- 
tially weakly interacting with the chaotic NMs in the N 
part, their excitation - if at all - will take times which 
increase exponentially with growing distance. Therefore, 
the wave packet will spread (if at all) slower than any 
power law. Thus, the LNL model is the only model 
we consider here, where almost trivial slowing down of 
spreading is expected. 

Recently, the dynamics of a similar system was theo- 
retically investigated in [24j , where a disordered subsys- 
tem of coupled anharmonic oscillators, linearly coupled 
to an infinite harmonic system, was considered. In this 
work the conditions which permit the persistence of the 
discrete breather of the isolated anharmonic system for 
small but nonvanishing couplings to the harmonic lattice 
were derived, and cases characterized by energy transfer 
to the harmonic system were also discussed. 

The NLN chain is expected to behave differently. As 
long as the wave packet is confined mainly to the L region, 
the dynamics is regular, and no spreading should occur. 
For large enough time, some part of the packet will leak 
out into the N regions. Therefore, finally spreading of 
the wave packet should occur. 



B. Numerical results 

1. Methods 

We consider compact DNLS wave packets at t = 
spanning a width L centered in the lattice, such that 
within L there is a constant initial norm of rij n = 1 and a 
random phase at each site , while outside the width L the 
norm density is zero. In the KG case, we excite each site 
in the width L with the same energy, E = H / L, i.e. initial 
momenta of pi = ±\ // 2E with randomly assigned signs. 

We use symplectic integration schemes of the SABA 
family of integrators [HI, l25l . [26j for the integration of 
equations ([2]) and (|4]). The particular symplectic scheme 
used for the DNLS model is described in the Appendix. 
The number of lattice sites N and the integration time 
step t varied between N = 1000 to N = 2000 and 



r = 0.01 to t = 0.1, in order to exclude finite size ef- 
fects in the wave packet evolution, and in order to reach 
long integration times up to 10 7 — 10 9 time units with 
feasible CPU times. In all our simulations, the relative 
energy and norm errors are kept smaller than 10~ 3 . For 
each parameter set we averaged our data over 1000 dif- 
ferent disorder realizations, unless otherwise stated, and 
denote this by (...). In particular, we compute mi and 
P, and we smooth (log 10 7712) and (log 10 P) with a lo- 
cally weighted regression algorithm 2JJ , and then apply 
a central finite difference to calculate the local derivatives 
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2. Weak disorder 

In this subsection we considerably extend the reports 
on the observation of weak chaos, strong chaos, the 
crossover between both, and the selftrapping regime in 
Ref. [H. In Fig. Q] we show results for the DNLS model 
with W = 4 and L = V = 21, for six different val- 
ues of the nonlinearity strength f3. The time evolution 
of (log 10 m 2 (i)) and (log 10 P(£)) is plotted in Figs. QJa) 
and [TJb), respectively. The evolution of the compact- 
ness index (£(i)) is shown in Fig. [He). In Figs. (Ud) 
and[TJe) we plot the time dependence of the numerically 
computed derivatives a m (t) and ap (t) ([9]) obtained from 
the smoothed curves of Figs. QJa) andQJb) respectively. 
Finally, in Fig. [TJf) the values of (SV(i)) are plotted. 

The weak chaos dynamics is observed in Fig. [1] for 
f3 = 0.012 (black curves) and /3 = 0.04 (magenta curves). 
Initially, the wave packets remain localized and all quan- 
tities of Fig. [T] are constant with a m , ap being practi- 
cally zero. After some detrapping time td the wave pack- 
ets start to subdiffusively spread with TO2 ~ t 1 ^ 3 and 
P ~ t 1 / 6 (Figs, [^a), (b), (d) and (e)). In addition, the 
compactness index (Q ss 3 (Fig. Q^c)), indicating that 
wave packets are well thermalized inside. The tendency 
towards complete derealization of wave packets is clearly 
depicted in the evolution of the averaged norm fraction 
(Sv) which remains at the L = 21 initially excited sites 
(Fig. [T](f )). After the detrapping time td (Sy) decreases 
continuously up to the final integration time t = 10 7 . 

Increasing the value to /3 = 0.18, the initial spreading 
dynamics enters the crossover between weak and strong 
chaos, and a faster spreading is observed. Spreading sets 
in earlier, the compactness index again indicates ther- 
malized wave packets, and the local derivatives a m and 
ap increase up to 0.4 and 0.2 respectively, with a pos- 
sibly very slow decreasing at even larger times. (5V) 
again continuously decreases to zero indicating complete 
derealization. 

For j3 = 0.72, we fully enter the strong chaos regime. 
Most importantly we observe a saturation of the local 
exponent a m around the theoretical value 1/2, with a 
subsequent deca y, ag ain as predicted by theory [20( , and 
first observed in 11911. 
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FIG. 1: (Color online) DNLS, W = 4: Evolution of (a) (log 10 m 2 (t)), (b) (log 10 P{t)}, (c) (C(<)), (d) finite-difference derivative 
a m (t) for the smoothed m2 data of panel (a), (e) finite-difference derivative ap(t) for the smoothed P data of panel (b), and 
(f) (Sv(t)) for the spreading of wave packets with initial width L = 21 and /3 = 0.012, 0.04, 0.18, 0.72, 3.6, 8.4 [(bl) black; 
(m) magenta; (r) red; (b) blue; (g) green; (br) brown]. In panels (a), (b), (d) and (e) straight lines correspond to theoretically 
predicted power laws m.2 ~ t a , P ~ t a ^ 2 with a = 1/3 (dashed lines) and a = 1/2 (dotted lines). Error bars in panels (a) and 
(b) denote representative standard deviation errors. 



Finally, for j3 — 3.6 and /3 = 8.4 (green and brown 
curves respectively in Fig. QJ the dynamics enters the 
selftrapping regime. We observe that a part of a wave 
packet remains localized, while the remainder spreads. 
The spreading portion results in a continuous increase of 
?7i2 (Fig. HJa)) which initially is characterized by large 
values of a m > 1/2 (Fig. Hfd)). For larger time a m de- 
creases below 1/2. The evolution appears to be rather 
complex, and is not captured by the theoretical consider- 
ations in [2(J ■ The large values of a m may be due to tem- 
poral trapping and detrapping processes in this strongly 
nonequilibrium dynamics of the wave packet, with a fi- 
nal trapped packet fraction remaining. Therefore, we 
observe that (log 10 P) starts to level off (Fig. QJb)), ap 
tends to very small values (Fig. [He)) and (£) tends to 
zero (Fig. (He)). In Fig. [ljf), we observe that the values 
of (Sy) saturate to higher values for f3 = 8.4. 

For these typical cases of weak chaos (f3 = 0.04), strong 
chaos (f3 = 0.72) and selftrapping ((3 = 3.6) we present 
in Figs. |5fa)-(c) the time evolution of the averaged norm 
density distributions (zi) in real space. All simulations 
presented in Fig. [2] started from the same initial profile 
with size L = V, therefore the width of the localization 
volume set by the linear case is the width of the distri- 
butions at the shortest times in the plots. In the weak 
chaos regime (Fig. [UJa)), the wave packets remain close 
to their initial configuration for some times, as the high 
(zi) values in the region of the initial excitation indi- 



cate, followed by derealization at larger times. Thus, at 
t = 10 7 the averaged wave packet has spread to about 
300 sites with (z t ) > 10~ 5 . Therefore, the wave packet 
spreads continuously over distances which are an order 
of magnitude larger than the limits set by the linear the- 
ory and destroying Anderson localization. In the case of 
strong chaos (Fig.f2Jb)), spreading is even faster, leading 
to more extended profiles at t = 10 7 : about 700 sites 
with (zi) > 10" 5 . In the selftrapping regime (Fig. f^c)), 
the spreading part of the wave packet covers 1500 sites, 
with another clearly visible part staying selftrapped at 
the initial excitation region. The curved fronts in the 
density plots in Fig. [2] follow from the theoretical pre- 
diction mi ~ t a with which leads to a packet width 
A/" - y/rn2 and Af ~ c M°gio *)/2. 

For the KG model ([3]) with W = 4 we present in 
Fig. [3] similar results to the ones for the DNLS model. 
For small values of the initial energy density E = 0.01, 
the characteristics of the weak chaos regime are observed: 
?7i2 ~ t 1 / 3 (black curve in Fig. |3{a)) after a detrapping 
time td ~ 10 5 , wave packets remain compact as they 
spread since (£) ps 3 (Fig.[3[b)), and the fraction (Hy) of 
the energy of the initially L = 21 excited sites decreases 
(Fig. Eljd)). For E = 0.04 we enter the crossover region 
between the weak and the strong chaos regimes, with 
characteristics similar to the DNLS case. For E = 0.2 
(red curves in Fig. 13]) we observe the typical behavior of 
the strong chaos scenario: spreading is characterized by 
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a saturated a m ~ 1/2 (Fig. [3jc)) for about two decades 
(log 10 1 « 3.5 — 5.5), followed by a crossover to the weak 
chaos dynamics with a m decreasing. Getting closer to 
the selftrapping regime for E = 0.75 (blue curves in 
Fig. [3]), or being deep inside it for E = 3 (green curves 
in Fig. [3]) the characteristics of the selftrapping behavior 
appear, since (£) decreases (Fig. [3jb)), and (Hy) tends 
to stabilize to non-zero small values (Fig. Efd)). Similar 
to the (3 = 3.6 and f3 = 8.4 cases of the DNLS model, the 
evolution of m-i (Fig. |3Ja)) shows an initial phase of fast 
growth with a m > 1/2 (Fig.|3jc)) followed by a lowering 
in the values of a m . 

Our numerical results are in accord with the predic- 
tions of weak and strong chaos regimes, as well as of 
the crossover from strong to weak chaos. Selftrapping is 
observed as well, with less understood strongly nonequi- 
librium dynamics of the trapped and spreading packet 
parts. We vehemently stress that in all our simula- 
tions we never observed any evidence of a wave packet 
transition from the weak chaos regime, characterized by 
a m = 1/3, to a subsequent slowing down of spreading, 
which would lead to a m < 1/3. 



3. Strong disorder 

To search for potential deviations from the predicted 
spreading laws, we turn to large values of W for the 
DNLS system. In all our simulations we had L > V . 
In particular, we set L = 10 and considered the cases 
with W = 15 and W = 40. For large values of W we ex- 
pect only the weak chaos and the selftrapping dynamical 
regimes to be observed (T^|. For each value of W three 
different values of /3 were considered, one being in the 
weak chaos regime and the other two in the selftrapping 




i°g w t log 10 t 



FIG. 3: (Color online) KG, W = 4: Evolution of (a) 
(Iog 10 m a (t)>, (b) (C(t)>, (c) a m (t), and (d) (H v (t)) versus 
log 10 t for the spreading of initially compact wave packets of 
width L — 21 with E = 0.01, 0.04, 0.2, 0.75, 3 [(bl) black; (m) 
magenta; (r) red; (b) blue; (g) green]. In panels (a), and (c) 
straight lines correspond to the theoretically predicted power 
laws 77i2 ~ t a with a = 1/3 (dashed lines) and a — 1/2 
(dotted lines). Error bars in panel (a) denote representative 
standard deviation errors. 



regime. In particular, we considered j3 = 0.5, (3 = 9 and 
ji = 30 for W = 15 and P = 1, = 25 and (3 = 100 for 
W = 40. The obtained results are shown in Fig. [4] for 
W = 15 and in Fig. [U for W = 40. 

For large W, the localization volume V decreases dras- 
tically, so that eventually the overlap integrals become 
small as well. Therefore, the values of /3 at which spread- 
ing becomes visible will increase. In the weak chaos 
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FIG. 4: (Color online) DNLS, W = 15: Evolution of (a) 
(log 10 m 2 (t)), (b) (£(t)), (c) a m (t), and (d) (S L (t)) versus 
log 10 £ for the spreading of initially compact wave packets of 
width L — 10 with /3 = 0.5, 9, 30 [(m) magenta; (g) green; 
(br) brown]. Mean values are averaged quantities over 100 
disorder realizations. In panels (a), and (c) straight lines 
correspond to theoretically predicted weak chaos behavior 
m 2 ~ t 1/3 . 
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FIG. 5: (Color online) DNLS, W = 40: Evolution of (a) 
(log 10 ma(*)>, (b) (at)), (c) am(t), and (d) (,&(*)) versus 
log 10 t for the spreading of initially compact wave packets of 
width L — 10 with f3 = 1, 25, 100 [(m) magenta; (g) green; 
(br) brown]. Mean values are averaged quantities over 100 
disorder realizations. In panels (a), and (c) straight lines 
correspond to theoretically predicted weak chaos behavior 
m 2 ~ t 1/3 . 



spreading regime the detrapping time t d becomes large 
and we start observing spreading after long time inter- 
vals. For W = 15 and /3 = 0.5 we find t d 10 3 
(Fig. EJa)). The local derivative a m increases from zero 
showing a tendency to approach the theoretically pre- 
dicted value a = 1/3 (Fig.|t[c)), and both (Q (Fig.^b)) 



and (Sl) (Fig.[4jd)) start to decrease. We note that since 
V ~ 1 for large W, we measure the time evolution of the 
fraction Sx(£) of the norm density of the L = 10 initially 
excited sites. 

The detrapping time td increases as W increases. This 
is seen from the results for W = 40, (5 = 1 (magenta 
curves in Fig. [5]). In this case, we have to wait at 
least up to td = 10 5 in order to get some evidence that 
spreading starts, since after that time (log 10 7712) starts to 
slightly grow (Fig. Efa)), while (C) (Fig.^b)) and (S L ) 
(Fig. EJd)) start to decrease. This increase of td happens 
despite the fact that the nonlinearity strength /3 also in- 
creased by a factor of two as compared to the W = 15 
case. Nevertheless even in this case of large W = 40 we 
are able to numerically observe the onset of spreading in 
the weak chaos regime. Increasing W to even higher val- 
ues pushes td to values larger than the final integration 
time t = 10 7 used in our simulations. 

With increasing (3 we observe selftrapping, but a part 
of the wave packet spreads and mi increases (green and 
brown curves in Figs. HJa) and[5]Ja)). The average com- 
pactness index (£) decreases, which is a clear indication 
that a part of the wave packet remains localized, and 
reaches smaller final values for larger j3 (Figs. 0|b) and 
(5Jb)). The selftrapping of the wave packets is also clearly 
seen from the evolution of (Sl}- In Fig.[4jd) we see that 
for W = 15, j3 = 9 (green curve) and /3 = 30 (brown 
curve) (Sl) decreases due to the spreading of a part of 
wave packets, while, finally it shows a tendency to level 
off to a positive value, indicating that part of the wave 
packets remains localized. Similar behaviors of (Sl) are 
observed in Fig. E^d) for the W = 40 case with j3 = 25 
(green curve) and (3 = 100 (brown curve), although the 
plateauing of (Sl) is not as clear as in Fig.[4jd). The nu- 
merically computed exponents a m exhibit the typical be- 
havior of selftrapping seen in Fig.QJd). For W = 15 they 
increase reaching values larger than 1/3 and afterwards 
decrease towards a m = 1/3 (Fig. Hfc)). F° r W — 40 a 
similar behavior is observed for f3 = 100, while for /3 = 25 
a m seems to approach the theoretically predicted value 
1/3 from below. 

Therefore, even for strong disorder, the dynamics of 
wave packets evolves according to the theoretical pre- 
dictions. Most importantly, we do not observe a slowing 
down of the wave packet below the limits set by the weak 
chaos regime. 



4. FSW model 

To further probe a possible slowing down of wave 
packet spreading beyond the limits of weak chaos, we turn 
to the FSW model ©. In Fig. [HI we present results with 
initial compact wave packets of width L = 21, with en- 
ergy density E = 0.05, similarly to the KG model ([3]). Wc 
observe that also for this model subdiffusivc spreading oc- 
curs, because the second moment and the participation 
number (red and blue curves respectively in Fig. Eta)) 
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FIG. 6: (Color online) FSW: Evolution of (a) (log 10 m 2 (t)> 
[(r) red curve] and <log 10 P(t)) [(b) blue curve], (b) (((t)), 
(c) a m (t) [(r) red curve] and otp{i) [(b) blue curve], and (f) 
(HL(t)) versus log 10 t for initially compact wave packets of 
width L — 21 with E = 0.05. 

increase continuously, and the corresponding exponents 
a m and ap (Fig. Etc)) tend to eventual constant non- 
zero values. The fraction (Hi,) of energy remaining in 
the L — 21 initially excited sites (Fig. HJd)) decreases as 
time increases indicating the dclocalization of the wave 
packets. The compactness index (Fig. E^b)) has a dif- 
ferent behavior with respect to what we have seen in 
the rest of our simulations, as it decreases slowly up to 
t « 10 , with a subsequent increase. Therefore, the wave 
packet is highly inhomogeneous for almost all the inte- 
gration time, violating the assumptions which are used 
in the theoretical considerations of weak chaos [2(| . The 
observed subdiffusive spreading may still not be in its 
final asymptotic range. Still, we again do not see any 
signature of a slowing down of this subdiffusive process, 
as was reported in [28| where a similar model was con- 
sidered. Clearly the FSW model calls for a thorough and 
independent study. 



5. LNL and NLN models 

In Fig. [7] we present results for the KG, LNL and NLN 
models for W = 4, L = 21 and E = 0.02. For com- 
parison we also include the results for the KG model ^ 
(magenta curves) with E = 0.02, for which subdiffusive 
spreading in the weak chaos regime is observed. The 
KG-NLN model (green curves in Fig. E|) exhibits a sim- 
ilar behavior, since both the second moment (Fig. Eta)) 
and the participation number (Fig. Etb)) start to grow 
after some detrapping time td rj 10 5 . This time is larger 
than the detrapping time of the KG model (td rj 10 4 ), 
because a wave packet in the KG-NLN model initially 
evolves in an almost linear system and only after some 
large time, when it has spread significantly to the nonlin- 



ear part of the lattice, spreading takes on characteristics 
of the purely nonlinear model. 

On the other hand the evolution of all quantities of 
Fig. El for the KG-LNL system (red curves) follows the 
KG model until t« 10 4 , because initially the wave pack- 
ets evolve in the same nonlinear system. Later on the 
wave packet enters the L (linear) parts of the system. 
Thus, spreading starts to retard, and both (log 10 777,2(7;)) 
(Fig. Eta)) and (log 10 P(t)) (Fig.[7Jb)) show a character- 
istic slowing down in the exponents a m (Fig. Et°0) and 
ap (Fig. Et e ))- In addition, (Sy) (Fig. Elf)) saturates at 
finite non-zero values, indicating that wave packets tend 
to localize again. For all three KG models, the values 
of (() (Fig. Efc)) show that wave packets do not become 
sparse and inhomogeneous in the course of time. We ob- 
tained similar results for the LNL and NLN DNLS models 
for W = 4, L = 21 and /3 = 0.04. 

Thus, also for the LNL and NLN models spreading is 
observed, and only in the case of the LNL models we have 
observed a slowing down of the spreading, as expected. 



IV. SUMMARY AND CONCLUSIONS 

We considered several models of disordered nonlinear 
one-dimensional lattices and performed extensive numer- 
ical simulations of norm (energy) propagations. Since we 
focused on the dynamical spreading of fronts, we pre- 
pared initial block wave packet profiles, having widths 
equal or larger than the average localization volume de- 
fined by the linear problem. While not performed, we ex- 
pect similar behaviors for initial Gaussian profiles, again 
where the width (for Gaussians, say the standard devia- 
tion) is on par with the average localization volume. 

We carefully studied statistical properties of the dy- 
namics, by varying the values of disorder and nonlincarity 
strengths over a wide interval, and by averaging results 
over many disorder realizations. Our results agree quite 
well with our theoretical expectations for the existence 
of the weak and strong chaos regimes. 

The main outcome of our study is that in the presence 
of nonlincaritics we always observe subdiffusive spread- 
ing, so that the second moment grows initially as rri2 ~ t a 
with a < 1, showing signs of a crossover to the asymp- 
totic ?7i2 ~ t 1 ^ 3 law at larger times. Remarkably, sub- 
diffusive spreading is also observed for large disorder 
strengths, when the localization volume (which defines 
the number of interacting partner modes) tends to one. 
Frohlich-Spencer- Wayne models which take the disorder 
strength to its infinite limits, are also showing subdif- 
fusive growth. Most remarkably, in none of our studies 
(except the artificial LNL case) did we encounter a slow- 
ing down of spreading beyond the limits set by the weak 
chaos predictions. Therefore, our numerical data sup- 
port the conjecture, that the wave packets, once they 
spread, will do so up to infinite times in a subdiffusive 
way, bypassing Anderson localization of the linear wave 
equations. 
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FIG. 7: (Color online) Evolution of (a) (log 10 m 2 (t)), (b) (log 10 P(t)>, (c) (C(t)> , (d) a m (t), (e) a P (t), and (f) (H v (t)) versus 
log 10 i for the spreading of initially compact wave packets of width L = 21 with W = 4 and _B = 0.02 in the KG [(m) magenta 
curves], the KG-NLN [(g) green curves] and the KG-LNL [(r) red curves] models. In panels (a), (b), (d) and (e) straight lines 
correspond to theoretically predicted power laws mi ~ t 1 ^ 3 , P ~ t 1 ^ 6 of the weak chaos regime. Error bars in panels (a) and 
(b) denote representative one standard deviation errors. 



The only cases where spreading shows a tendency to 
stop are the LNL models, for which nonlincarities are ab- 
sent everywhere except inside a finite-size central region, 
where the initial wave packet is launched. In these mod- 
els, when wave packets have spread substantially, their 
chaotic component in the central region of the lattice 
becomes weak, and distant normal modes in the linear 
parts of the system are exponentially weakly coupled to 
the central nonlinear region. 

When the nonlinearity strength tends to smaller val- 
ues, waiting (detrapping) times for wave packet spreading 
of compact initial excitations increase beyond the detec- 
tion capabilities of our computational tools. The cor- 
responding question of whether a KAM regime can be 
entered at finite nonlinearity strength was addressed in 
[l2j and is analyzed in detail in a forthcoming work [29| . 
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Appendix A: Symplectic integration of the DNLS 
equations 

We first discuss a novel PQ method which we designed 
to integrate the DNLS equations locally. Previously used 
methods employ a transformation of the wave function 
from real into Fourier space and back, at each integration 
step. These transformations induce small but observable 
corrections in the tails of the wave packet, which slowly 
but steadily grow in time. In such a case we will have to 
stop the integration once this noisy background reaches a 
substantial level. The PQ method avoids the generation 
of this background by simply not performing the Fourier 
transformation. Instead the PQ method integrates the 
DNLS equations in real space. 

The canonical transformation 

in = ^{qi + ipi) (Ai) 
of the complex variable ipi in Eq. (Q} transforms (p} into 

= ^2^(qf +pf) + g(«f +pf) 2 - (qi+iqi+Pi+iPi). 

(A2) 

where qi and pi are generalized coordinates and momenta, 
respectively. 

If a Hamiltonian function can be split into two inte- 
grable parts, then a symplectic integration scheme can 
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be used for the integration of its equations of motion. 
One possible splitting of the DNLS Hamiltonian (|A2|) 
into two separate Hamiltonian functions A and B is 



.4 
D 



e l i 2 i 2-i , P f 2 , 2\2 
I +Pl)+ e(<ll +Pl) > 



^2 (9 ' ' 8 

i 

~22(Ql+lQl +P1+1P1) 



(A3) 



are intcgrablc, under the corresponding operators 



and 



„Tif 



Pi 
'i,' 



Pi 



(pj_i +p;+i)r 



(A6) 



Hamiltonian A is integrable and the operator e rL - 4 which 
propagates the set of initial conditions (qi,pi) at time t, 
to their final values (?» 3 pj) at time i + r is 



= qicos(aiT) 
= pi cos(air) 



-pi sin(a;r) 
- qi sin(azr) 



(A4) 



with a t = ei+P(qf+pj)/2. Hamiltonian B of Eq. (|A3|) is 
not integrable, thus the operator e rLB cannot be written 
explicitly. If we consider B as a separate Hamiltonian 
function and again split it as B = P + Q, the component 
parts 



p = -'Y^pi+xpu q = -^2 



qi+rn- 



(A5) 



Vi = 1i 

p'i = pi + (m-i + qi+i) T 



(A7) 



This technique of splitting the Hamiltonian into multiple 
parts has been used in different applications of symplectic 
integrators (see for example (30l|). 

In our simulations we successively apply the SBAB2 
symplectic integrator [15, 25, 26] twice: first for the split 
(Hd = A + B) of the DNLS Hamiltonian, Eq. (jAlJ) . and 
second for the split B = P + Q in Eq. (|A3]) . The solu- 
tion for the equations of motion from the Hamiltonian 
Eq. (|A2[) is thus approximated by the application of 13 
simple operators on an initial condition (qi,pi), since 



e d 1 TL Ae dic 2 TL P ^tLq e c 2 d 2 TL P ^tLq e dic 2 TL P e d 2 TLA e dic 2 TLp ^tLq e c 2 d 2 TLp e clrLq e d 1 c 2 TL P e diTL A 



with the SBAB 2 coefficients [23 of di = 1/6, d 2 = 2/3 and c 2 = 1/2. 
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